;;=====================================================================================
;; calculate RH pdf 
;;=====================================================================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "pdfx_log.ncl"
load "pdfxy_xlog.ncl"

begin

 klev = 20
 nbin = 25

 nbinx = nbin
 nbiny = nbin

 print("")
 print("")

 rcase = getenv("rcase")

 path = "./"

 print("")
 print("")

 fba = "icnc3h_sp_T63L31K01_2000_jan-jun.nc"
 fbb = "wsub3h_sp_T63L31K01_2010_jan-jun.nc"
 fbc = "st_T63L31K01_2010_jan-jun.nc"

 fna = path+fba
 fnb = path+fbb
 fnc = path+fbc
 print("input file :   "+fna)
 print("input file :   "+fnb)
 print("input file :   "+fnc)
 print("")
 print("")

 fla = addfile(fna,"r")
 flb = addfile(fnb,"r")
 flc = addfile(fnc,"r")


 xt   = flc->st
 xi   = fla->ICNC
 xh   = flb->W_TURB

 xi   = xi / 1000. ;; /m3 -> /L 

 xh   = where(xh.le.0.001,-999,xh)
;xh   = where(xh.gt.1.2,-999,xh)
 xh@_FillValue = -999

 ;; cm/s 

 xh   = xh * 100.

 ;; exclude mixed phase 

 xi   = where(xt.lt.238.,xi,-999)
 xi@_FillValue = -999


 print("read data :   done !")
 print("")

 ;;----------------------------------------------------------
 ;; reorder data 
 ;;----------------------------------------------------------

 sxt = ndtooned(xt(:,0:klev,:,:))
 sxi = ndtooned(xi(:,0:klev,:,:))
 sxh = ndtooned(xh(:,0:klev,:,:))
 ;;sxr = ndtooned(xr(:,0:17,:,:))
 ;;sfr = ndtooned(fr(:,0:17,:,:))


 sxh@_FillValue = -999
 sxt@_FillValue = -999
 sxi@_FillValue = -999

 print("60%")


 px = sxi  ;;(::10) 
 pa = sxh  ;;(::10) 

 print("60%")


;*************************************************************
; Bivariate (Joint) PDFs
; Use coordinate subscripting to select data
;*************************************************************

 opa         = True
 opa@bin_min =    0.
 opa@bin_max =   24.

 pdf1_0 = pdfx_log(sxh,nbinx,opa)


 opt          = True

 opt@binx_min =    0.1 ;;- minimum value for the x bin boundary.
 opt@binx_max =  10000 ;;- maximum value for the x bin boundary.
 opt@biny_min =   0 ;;- minimum value for the y bin boundary.
 opt@biny_max =  120 ;;- maximum value for the y bin boundary. 

 pdf2_0 = pdfxy_xlog(px,pa,nbinx,nbiny,opt)

 fno = "mod_sparticus_wsub_ni_pdf_25.nc"

 system("rm "+fno)

 flo = addfile(fno,"c")

 flo->pdf1_0=pdf1_0

 flo->pdf2_0=pdf2_0

end



